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Abstract 

We investigated numerically the high-temperature/high-pressure phase di- 
agram of Xenon as modelled through the exp-6 interaction potential, which 
is thought to provide a reliable description of the thermal behaviour of rare 
gases under extreme conditions. We performed a series of extensive NVT 
Monte Carlo simulations which, in conjunction with exact computation of 
the solid free energy by the Frenkel-Ladd method, allowed us to precisely lo- 
cate the freezing and the melting thresholds at each temperature. We find 
that, under isothermal compression, the exp-6 fluid freezes directly into a FCC 
solid; however, above 4500 K, an intermediate BCC phase becomes stable in 
a narrow range of pressures. The chemical potential of the HCP phase does 
never significantly differ from that of the FCC solid of equal T and P, though 
the former is found to slightly overcome the latter. We discuss our results in 
the light of previous numerical studies of the same model system and of the 
experimental data available for Xenon. 
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I. INTRODUCTION 



In the last decade, there has been an increasing interest in the high-temperature/high- 
pressure (HT/HP) properties of many materials [1-3], as realized for instance in the deep 
core of the planets of our Solar system. Such peculiar conditions do in fact provide a new 
means for testing current condensed-matter theories, since squeezing matter to planetary 
pressure and heating it up to a few thousand degrees can trigger various forms of structural 
reorganization both at the macroscopic and at the molecular level (see, for instance, the 
two paradigmatic cases of water and methane), or even cause modifications in the electronic 
transport properties (as it happens for e.g. hydrogen). Since down-the-earth experiments 
are somehow limited (pressure can hardly be pushed over a certain threshold and, more 
important, any huge compression is plagued by severe non-hydrostaticity problems), the 
only possible insights into the transformations undergone by many substances at extreme 
conditions do often come from numerical simulation - provided, however, a theoretical model 
of that substance is amenable to careful investigation and numerical errors are kept under 
control. 

Due to their closed electronic shells, (heavy) rare gases are usually believed to be the 
simplest substances of all, hence it may appear somewhat strange that their thermal be- 
haviour in the HT/HP regime (i.e., in the very-dense fluid and solid regions) is not yet well 
assessed [4]. Despite the fact that one can get rid of quantum-mechanical considerations 
almost completely (at a temperature as high as 3500 K, it was shown in Ref. [5] that elec- 
tronic excitations do not play any significant role up to pressures of order 50 GPa), a classical 
interaction potential that accurately reproduces the thermodynamic properties of rare gases 
cannot be of the simple Lennard- Jones form. In fact, while performing very well at ambient 
conditions, the Lennard- Jones potential loses much of its reliability when temperature and 
pressure take huge values. For high system densities, three-body contributions to the effec- 
tive potential are likely to be as important as the two-body term [6], with the effect that 
the atomic core is softer than implied by the Lennard- Jones interaction. If one insists in 
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using a two-body effective potential, a more reliable form for rare gases turns out to be the 
(modified) Buckingham or exp-6 potential [7] (see it defined at Eq. (2.1)), as parametrized 
through high-density experimental data [8]. As a matter of fact, when HT/HP conditions 
hold, an exponential law is a more adequate representation for the interatomic repulsion at 
short distances than is a power law. In this respect, we can say that the exp-6 potential 
takes into account the effects of the three-body interaction in an isotropic fashion. 

Recently, some controversy has arisen about the topology of the HT/HP part of the phase 
diagram of Xe, whose behaviour should be representative also of Ar and Kr. Stimulated 
by the findings of a recent laser-heated diamond-anvil-cell (DAC) experiment on Xe [9], 
Belonoshko et al. have done molecular-dynamics (MD) computer simulations of the exp-6 
potential in order to interpret those data [5,10]. In the DAC experiment, a Xe fluid sample 
is compressed at a given pressure in the range between 10 and 40 GPa and subsequently 
heated with a laser until a bulk transition is detected. The experiment is meant to provide 
information about the freezing of fluid Xe in the HT/HP region and, ultimately, to draw 
the fluid-solid coexistence line in the P-T diagram. Surprisingly, the freezing line shows an 
unexpected cusp near T = 2700 K and P = 15 GPa, which was originally imputed to the 
boundary between a fluid-solid and a fluid-glass transition [9]. The numerical simulations of 
the exp-6 model by Belonoshko et al. have instead revealed a competition between two dis- 
tinct, FCC and BCC solid phases, which happen to exchange their relative thermodynamic 
stability right in the region of the observed cusp. More precisely, for temperatures above 
T ~ 2700 K, a fluid-BCC-FCC sequence of phases is reported upon isothermically increasing 
the pressure beyond 25 GPa. This interpretation has been criticized by Kechin [11]. 

The simulation method used in Refs. [5,10] was the so-called two-phase or coexistence 
method [12,13], as implemented into an isothermal-isobaric MD simulation [14]. In point 
of principle, this technique should be able to detect a coexistence between two structurally 
distinct phases whenever it occurs; in practice, this method has some limitations since it 
requires large sizes and very long simulation times to reach good equilibration in the region 
of the interface between the coexisting phases [15]. Moreover, it cannot be excluded that the 
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route to equilibrium of an unstable interface will pass through (and long remain in) a phase 
which, though not being the preferred one under bulk conditions, is nonetheless promoted 
by (small) spatial inhomogeneities of the system. This would be especially so when the 
difference in Gibbs free energy between two distinct solid phases is very small: in that case, 
a system being prepared in a metastable solid phase would take very long time (possibly 
much longer than accessible to simulation) to relax to the truly stable solid. 

In order to gain further insight into the nature of Xe freezing at the temperatures in- 
vestigated by the DAC experiment, we have carried out a Monte Carlo simulation in the 
canonical ensemble of the same exp-6 system of Refs. [5,10], but now calculating the free 
energy of the relevant solid phases by the method of Frenkel and Ladd [16], which would 
provide the benchmark for solid- free-energy evaluations. When used in conjunction with 
conventional thermodynamic integration, this method would allow one to draw the "exact" 
phase diagram of the given model system and the only source of error would be associated 
with the finite size of the sample. Indeed, if the chemical-potential gaps between different 
solid structures are minute, also the statistical errors affecting the relevant thermodynamic 
averages are to be made sufficiently small (i.e., by carrying out long enough simulation 
runs). As a matter of fact, we do neither confirm the findings of Belonoshko et al. - for 
we predict that the BCC solid would become stable in Xe only above 4500 K - nor find 
any good agreement with the DAC data of Ref. [9] at the highest temperatures, i.e., no 
clear bump is seen along the freezing line. A similar outcome was also found by Frenkel 
for He [17], whose BCC-FCC transition would only occur in a region of the phase diagram 
where the fluid is stable. 

The outline of the paper is the following: in Section 2, we introduce the model system 
and describe our simulation method. In Section 3, we present our numerical results and 
draw the ensuing phase diagram. Then, in Section 4, we make a critical comparison of our 
outcome with those of previous studies. Further comments and remarks are postponed to 
the Conclusions. 
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II. MODEL AND METHOD 



A. The exp-6 model 



All textbooks in statistical mechanics report that rare-gas thermal behaviour is well 
accounted for by the simple Lennard- Jones pair potential. However, it is less known that, 
when rare gases are very dense, a different potential with a softer repulsive shoulder is better 
suited. Such is the exp-6 potential, defined to be 



where a controls the softness of the repulsion and e > is the depth of the potential 
minimum located at r m . We select the value of a in such a way that, for the given a and 
r m , the function appearing in the second line of Eq. (2.1) reaches its maximum right at a. 
Hence, as r moves down to a, v (r) reaches a stationary value and then goes abruptly to 
infinity. Ross and Mc Mahan have shown that, for a suitable choice of the parameters a, e, 
and r m , v(r) gives a rather faithful representation of many thermodynamic properties of Ar 
at high densities. Moreover, a corresponding-state theory appears to hold, which permits 
to derive the model parameters appropriate to Kr and Xe from those of Ar. For Xe, such 
best parameters are a = 13, e/ks = 235 K, and r m = 4.47A [8], where ks is the Boltzmann 
constant (for a = 13, the value of a is then 0.246972 ... r m ). In the following, we work 
with dimensionless variables, i.e., T* = k B T/e, P* = Pr^/e, and p* = pr z m (p = N/V 
is the particle- number density). The accuracy of the exp-6 modelization for Xenon at low 
pressures was tested by Belonoshko et al. against available experimental results and, up to 
0.7 GPa, was found to be good (see Fig. 7 of Ref. [5]). 

Recently, the vapour-liquid equilibria of the exp-6 fluid were studied numerically through 
a series of Monte Carlo (MC) simulations in the Gibbs ensemble [18]. As a result, the critical 
point was located at p* = 0.303 and T* = 1.316. From a simulation study by Errington and 
Panagiotopoulos a slightly different estimate of the critical-point coordinates is extracted, 




(2.1) 
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namely p* = 0.318 and T* = 1.318 [19]. Finally, the freezing density of the exp-6 fluid was 
estimated for a number of temperatures by the Hansen- Verlet criterion [20,21]. 

B. Details of the Monte Carlo simulation 

The numerical-simulation method provides virtually exact information on the statistical 
behaviour of a given model system. We have performed Monte Carlo (MC) simulations 
of the exp-6 model in the canonical ensemble (i.e., at constant temperature T, volume 
V, and number N of particles), using the standard Metropolis algorithm for sampling the 
equilibrium distribution in configurational space. The values of N that we only consider are 
those that fit to a cubic simulation box with an integer number of cells, i.e., N = 4n 3 for 
the FCC solid and N = 2n 3 for the BCC solid, n being the number of cells along any spatial 
direction. For a given particle number, the length L of the box is adjusted to a chosen 
density value p, i.e., L = (N/p) 1 ^ 3 . Called a the distance between two nearest-neighbour 
reference lattice sites, we have a = (y/2/2) (L/n) for a FCC crystal and a = (v / 3/2)(L/n) 
for a BCC crystal. Finally, periodic conditions are applied to the box boundaries. 

Our largest sample sizes were N = 1372 (FCC) and iV = 1458 (BCC) which, to all 
practical purposes, can be regarded as nearly the same size. We point out that comparing 
the statistical properties of similar (large) sizes is mandatory when we want to decide which 
solid phase is stable at given T and P, since this produces comparable statistical errors for 
the two phases. Otherwise, one can resort to some extrapolation to the thermodynamic 
limit which is however computationally far more demanding. 

Some X-ray diffraction studies indicate that the FCC-HCP transition could be a common 
behaviour in all rare-gas solids [22-24]. At room temperature, the HCP solid would be stable 
above 75 GPa, while metallization is not expected before « 130 GPa [25]. We have thus 
carried out also a number of simulations of a HCP crystal hosting 10 x 12 x 12 = 1440 
particles (the simulation box is only approximately cubic in this case and the z axis is 
now oriented along the [111] direction - the three box sides have lengths of L x = n x a, 
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Ly = (y/3/2)n y a, and L z = (y/E/3)n z a, with n x = 10, n y = n z = 12, and a = (y^/p) 1 / 3 ). 

To locate the transition point at a given T, we have numerically generated two isother- 
mal quasistatic paths, starting from the very dilute fluid on one side and from the highly 
compressed solid on the other side. As a rule, in the solid region, the last MC configuration 
produced at a given p serves, after suitable rescaling of particle coordinates, as the starting 
configuration for the run at a slightly lower density. Similarly, in the dense-fluid region, the 
simulations are carried out in a chain, i.e., the run at a given density is started from the 
last (rescaled) configuration produced at a lower p value. The FCC and BCC solid paths 
are followed until the fluid spontaneously forms during the MC run, as evidenced by the 
abrupt change in energy and pressure. Usually, we were able to overheat the solid for a 
little beyond the fluid-solid phase boundary, while undercooling of the fluid is much easier. 
For each p and T, equilibration of the sample typically takes 2 x 10 3 MC sweeps, a sweep 
consisting of one attempt to sequentially change the position of all particles. The maximum 
random displacement of a particle in a trial MC move is adjusted once a sweep during the 
run so as to keep the acceptance ratio of the moves as close to 50% as possible, with only 
small excursions around this value. 

For given NVT conditions, the relevant thermodynamic averages are computed over a 
trajectory whose length ranges from 2 x 10 4 to 6 x 10 4 sweeps. The excess energy per particle 
u ex , the pressure P, and (in the solid phase) the mean square deviation 5R 2 of a particle 
from its reference lattice position are especially monitored. Pressure comes from the virial 
formula, 

P = pk B T+ { -^, V = -Wn,v\r lj ) (2.2) 

6 i<3 

(rij is the distance between particles i and j). In practice, to avoid double counting of 
interactions, the pair potential is truncated above a certain cutoff distance r c , which is 
taken to be only slightly smaller than L/2. Then, the appropriate long-range corrections 
are applied to energy and pressure by assuming g(r) = 1 beyond r c , g(r) being the radial 
distribution function (RDF). 
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The RDF histogram is constructed with a spatial resolution of Ar = r m /50 and updated 
every 10 MC sweeps. The RDF is computed up to a distance of L/2: at that distance, the 
g(r) was never found to significantly differ from unity, at least for the largest system sizes. 

To evaluate the numerical errors affecting the main statistical averages, we divide the 
MC trajectory into ten blocks and estimate the length of the error bars as being twice 
the empirical standard deviation of the block averages from the mean (under the implicit 
assumption that the decorrelation time of any relevant variable is less that the size of a 
block). Tipically, the relative errors of energy and pressure are of few tenths of a percent. 

The difference in excess free energy between two equilibrium states of the system, say 1 

and 2, lying within the same phase is computed through the combined use of the formulae 

/cx(T 2 , p) _ U(T uP ) fT2 ^ Ucx (T, p) 
T 2 Ti 7ti T 2 

and 

rp2 1 

PU(T, P2 ) = (3f cx (T,Pi)+ dp- 

Jpi p 

where f ex is the excess Helmholtz free energy per particle and (3 = (/c B T) _1 . The integrals in 
Eqs. (2.3) and (2.4) are performed numerically by applying the Simpson rule to a linear-spline 
approximant of the simulation data for energy and pressure. 

The above formulae are however useless if one does not have an independent estimate 
of the system free energy in a reference state. Only in this case, Eqs. (2.3) and (2.4) help 
finding the free energy of any other state in the same phase. The choice of such a reference 
state is different for the fluid and for the solid phase. As a reference state for the fluid, we 
can choose any equilibrium state been characterized by a very small p value and arbitrary 
T (say, a nearly ideal gas), since then the excess chemical potential of the system can be 
accurately estimated by the Widom or particle-insertion method, 

A* ex = ~k B T In (exp (-/3£ ins )) , (2.5) 

where Ei ns comprises all interaction terms between a randomly-inserted ghost particle and 
all the system particles. The average in Eq. (2.5) is evaluated numerically during a run of 
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(2.4) 



typically 5 x 10 4 equilibrium sweeps, with an insertion attempted at the completion of every 
sweep. Once /x ex is given, the excess values of free energy and entropy will follow from 

Pfcx = /9/iex - — + 1 and ^ = (3(u cx - / ex ) . (2.6) 

It is useful to note that, from a strictly numerical point of view, choosing a very dilute gas 
as a reference state for the fluid is far better than starting the thermodynamic integration 
in Eq. (2.4) from the ideal gas of equal temperature. In fact, unless one has a lot of thermo- 
dynamic points in the very-dilute region of the phase diagram, a spline interpolant of flP/p 
that is sufficiently accurate in this region is hard to construct. 

C. The Frenkel-Ladd calculation of solid free energies 

The Frenkel-Ladd (FL) method for calculating the free energy of a solid system relies on 
a different kind of thermodynamic integration [16,26]. The idea is to continuously transform 
the system of interest with potential energy V\ into a reference solid system of known free 
energy (typically, an Einstein crystal). This is accomplished through a linear interpolation 
of the potential-energy functions V and V\ of the two systems, i.e., V\ = V + A(Vi — Vo) 
(with < A < 1). Then, the difference in free energy between two homologous NVT states 
of the systems is calculated via the exact Kirkwood formula 

F!-F = f'dXiVi- V ) x , (2.7) 

J 

where the average (. . .) A is taken over the equilibrium distribution of the hybrid system with 
V\ potential and is evaluated numerically by a Monte Carlo procedure. 

Upon denoting {R^, i = 1, . . . , N} the reference lattice positions, the Einstein model is 
described by an interaction potential of 

1 N 2 

\/ Ein (R iV ) = -cE(R-Rf ) ) , (2.8) 
A i=i 

where c is the spring constant - the same for all oscillators. It readily follows from Eq. (2.8) 
that the free energy and mean square separation of a particle from its reference lattice site 
are given by 
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3N , f 2n \ , 3 



where A is the thermal wavelength. In deriving the first of Eqs. (2.9), no Gibbs factor of 
AM was included in the partition function since the Einstein oscillators are distinguishable 
entities. 

A non trivial problem with this choice of reference system, already pointed out in the 
original article [16], is the following: while the Einstein particles can only perform limited 
excursions around their reference positions, which implies a finite value of 5R 2 , there is no 
means to constrain particles interacting through V\ to remain confined in the neighborhood 
of their initial positions, even in the solid phase. In other words, at variance with V^in, V\ is 
a translationally-invariant potential, with the result that the integrand in Eq. (2.7) shows a 
divergence, in the thermodynamic limit, at A = 1. To overcome this problem, Monte Carlo 
simulations of V\ are usually performed under the constraint of a fixed center of mass. We 
do not repeat here the entire analysis showing how to deal with such a constraint in the 
equilibrium sampling, but simply quote the final result (see Ref. [26] for details): 



PU = P- 



, ( to s-zn\ 3. , n , , 2 In AT ln(2vr) In Upc)^ 2 ) 

+ ££d\(AV)™, (2.10) 

where the superscript CM denotes a constrained average and AV = V\ — V . Poison et al. 
have conjectured that, on fairly general grounds, (3f ex (N)+ln N/N = /3f eK (oo)+0(N~ 1 ) [26]. 
We shall see later whether this is found in our case. Given the value of / ex , the excess chemical 
potential of the solid still follows from the first of Eqs. (2.6). 

A few considerations on the numerical implementation of the FL method are now in order. 
While in principle the value of c can be chosen arbitrarily, just for numerical purposes it is 
convenient to take it such that SR 2 of the target solid is close to 3 /((3c): in this case, the 
variations of the integrand in Eq. (2.7) with A are likely to be small as well as the error in 
performing the numerical quadrature. 



A further caveat must be added for hard-core potentials (as is the one we are interested 
in here). If we write the potential as the sum of a hard-sphere part plus a remainder, 
V\ = Vhs + K, then it is mandatory to include a Vhs term also in Vb, since otherwise there 
is no way to transform continuously from V\ to Vq. As a result, 

V x = V HS + V Ein + A(K - V Ein ) , (2.11) 

implying AV — V r — Vein in Eq. (2.10). However, the free energy of a system of hard-core 
(that is, mutually interacting) Einstein particles is not known; what we can only say is 
that this interaction cannot have any consequence on the thermodynamic quantities of an 
Einstein crystal in so far as c takes sufficiently large values, since in this case particles are 
blind to each other. Then, unless c is given a very huge value, the FL calculation is only 
viable for a very cold solid. 

Whether a value of c is huge or not can be decided from the comparison between ^j3/(/3c) 
and (a — a)/2. Should the former be much smaller than the latter, the error committed in 
assuming the form (2.9) for the free energy of the interacting Einstein crystal is presum- 
ably negligible. Frenkel and Ladd have indicated a more rigorous method to ascertain the 
importance of hard-core interactions for the free energy of an Einstein solid, based on the 
estimate of the leading (virial) correction to the free energy of the non-interacting Einstein 
particles [16]. 



III. RESULTS 
A. Test calculations 

We have first checked our FL simulation code against the hard-sphere calculation given 
in Ref. [26]. We take a system of N = 256 hard spheres of diameter a at pa 3 = 1.0409 
(this is about the melting density of the FCC solid), with c = 6000 k B T/a 2 . We estimate 
numerically the average of AV = — V^in over the canonical distribution relative to V\ for a 
number of A values in the range from to 1. Precisely, we take A to increase in steps of 
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0.05 from to 0.9, with a smaller step of 0.01 in the range from 0.9 to 1. For each value 
of A, as many as 2 x 10 4 MC sweeps are produced at equilibrium. The MC moves are such 
that the position of the system center of mass is fixed; to the same purpose, if a particle 
happens to move (slightly) out of the simulation box, no attempt is made to put it back into 
the box [27]. Besides the value of (AV) X , we also estimate the statistical error affecting this 
quantity by partitioning the MC trajectory into large blocks. Finally, the integral over A in 
Eq. (2.10) is calculated using the Simpson rule. For the chosen c value, the virial correction 
to the free energy of the Einstein crystal due to hard-core interaction between particles is 
negligible (/3A/ ex w lO" 8 ). We thus find (3f cx + lnN/N = 5.891(5) for the FCC solid, which 
is fully consistent with the result obtained by Poison et al. (see Fig. 2 of Ref. [26]). 

As a further check of our numerical code, we have considered a model system of soft 
spheres interacting through the repulsive potential v(r) = e(a/r) 12 . We take a system 
of N = 256 particles with k B T/e = 1, per 3 = 1.1964, and c = 500 e/a 2 . We estimate 
numerically the average of Ay = V\ — (yf^ + Vein) over the canonical distribution relative 
to V\ for a number of A values in the range from to 1 (V^ is the total potential energy 
as calculated for particles located in the respective perfect-lattice positions). For such AV , 
the quantity /3Vf®/N must be added to the r.h.s. of Eq. (2.10) to obtain (3f ex . For each 
value of A, as many as 5 x 10 4 MC sweeps are produced at equilibrium. In the end, we find 
f3f ex + In N/N = 9.208(2) for the FCC solid, which again perfectly agrees with Fig. 1 of Ref. 
[26]. 

B. The exp-6 simulation 

Moving to the exp-6 model, we have calculated the excess free energy of both the FCC 
and the BCC solid for a number of (p*,T*) pairs, using for the integral in Eq. (2.10) the 
same specifications as given above for the test calculation on hard spheres. Though one 
single free-energy calculation by the FL method would be sufficient for estimating the free 
energy of the solid in any other state, we have often found more practical to repeat the FL 
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calculation rather than generating a long thermodynamic path to a very distant point in 
the phase diagram. In Table 1, we collect the calculated values of the excess free energy 
per particle in a few state points, for three different solid phases of the exp-6 model (and 
for just the largest sizes). The numerical precision reported for each /3/ ex is the Simpson 
integral of the errors that are associated with the values of (3 (AV)^ M /N, as calculated for 
the selected A points [28]. In the same table, we indicate the value of c* = cr^/e that is 
considered for each single state. In every case, this c is such that 3/ ((3c) ~ 5R 2 of the exp-6 
solid. Since all of the tabulated state points lie sufficiently far from the melting line, the 
average square excursion of an exp-6 particle from its reference lattice site turns out to be 
quite small. Hence, all of the c*'s are much larger than 1 and the virial correction to the 
Einstein-crystal free energy is absolutely negligible. 

We get an indirect check of the free-energy values in Table 1 by calculating, via ordinary 
thermodynamic integration based on Eqs. (2.3) and (2.4), the difference in excess free energy 
between various pairs of BCC and FCC states in the table. We have always found a perfect 
agreement with the tabulated values, to within the reported numerical precision, for both 
types of solids. For each (p*,T*) pair in Table 1, we have finally verified that the linear 
scaling of (3f ex (N) + In N/N as a function of N^ 1 holds well, as demonstrated in Fig. 1 in 
one case only. 

In order to see whether a stable BCC phase exists above T* ~ 11 (which is where the 
authors of Ref. [5] would locate the fluid-BCC-FCC triple point), we have carried out exten- 
sive MC simulations of the exp-6 model for a number of T* values (4.25, 8.15, 12.77, 16, 20, 
and 25). For all such values, we construct the fluid path and two solid, FCC and BCC paths 
(for T* = 4.25, only the FCC equation of state is generated). For T* = 16,20, and 25, we 
have also constructed the equation of state for N = 1440 particles in a HCP arrangement. 
With the only exception of T* = 4.25, where N = 864, the other isotherms are plotted for 
N = 1372 (fluid and FCC) and N = 1458 (BCC) particles. Once the free energy of the 
system is known along an isothermal path, the chemical potential fi along the same path 
readily follows from the first of Eqs. (2.6) as either a function of p or, equivalently, of P. 
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For given T and P values, the thermodynamically stable phase is the one with lower 
hence, the phase transition from fluid to e.g. FCC at constant temperature is located at the 
pressure where the chemical potential of the fluid takes over the \i of the FCC solid. 

Surprisingly, we find that, at variance with the conclusions of Ref. [5,10], the BCC phase 
of the exp-6 model is only stable at very high temperatures, i.e., above T* ~ 20, and in a 
narrow pressure interval. For all reduced temperatures up to 16, the scenario is similar (see 
it represented in Figs. 2 and 3 for T* = 16): as pressure grows, the chemical potential of 
the fluid eventually overcomes that of the FCC solid; at the crossing point Pfs, the latter is 
already lower than the BCC chemical potential. The BCC solid would be more stable than 
the FCC solid only for pressures lower than P ss (< Pfs), that is in a region where the fluid is 
the stable phase (however, at T* = 4.25, a pressure P S s simply does not exist, i.e., the BCC 
solid is never preferred to the FCC solid). Values of Pss, Pfs, an d of A/x = /ifcc ~ /^bcc ( as 
calculated at the fluid-solid transition pressure Pfs) are reported in Table 2 for the various 
temperatures. In the same table, the values of the freezing and the melting densities are 
also indicated. All the tabulated quantities are reported with three decimal digits, with 
no indication of the estimated statistical errors. The general question of the numerical 
reliability of our results is postponed to the next Section, but we can anticipate that the 
results are fully trustworthy. 

Fig. 2 shows the chemical potential excess A/x(P) of the FCC phase relatively to BCC, for 
T* = 16. This quantity is larger than zero (i.e., the stable solid would be of BCC structure) 
only inside the fluid region of the phase diagram. Hence, the fluid directly transforms into 
the FCC solid, while the BCC phase is metastable. In the same picture, we also plot the 
difference in chemical potential between FCC and HCP, a quantity which is nearly constant 
and close to —0.004 at all pressures. This means that the HCP solid is never stable, albeit 
it is about to be so at all pressures. 

Considering the small \x gaps between the various solids, one may wonder whether the 
above conclusions are influenced by the finite size of the simulated samples. To this purpose, 
we have taken T* = 16 and plotted in Fig. 2 the same quantity A/i(P) as before, but now 
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calculated for a smaller FCC solid of 864 particles and for a BCC solid of 1024 particles. 
The new curve differs from the old one just for a small rigid shift towards low pressures, and 
crosses zero at /3-Pss = 67.686 r~ 3 (to be contrasted with the previous value of 69.130r~ 3 ). A 
bit smaller is the variation of (3Pfs as N changes from 1372 to 864 (it moves from 71.877 r~ 3 
to 71.051 r m 3 ). We expect that similar results would be obtained at the other temperatures. 
Considering the large iV values involved, it is very unlikely that our conclusions for the 
largest sizes can be reversed in the thermodynamic limit. 

In Fig. 3, the various branches of the equation of state (3P vs. p are plotted for T* = 16. 
This curve shows a straight-line cut at P = Pps, which allows one to identify the values of pf 
and p m for the given T. Correspondingly, the Helmholtz free energy shows, as a function of 
the specific volume v — 1/p, a straight-line behaviour in the interval between p" 1 and pf 1 . 
This line is the common tangent of the fluid and the FCC free-energy curves (its equation 
is Pf = {3p>FS — fiPpsv)- Finally, a look at the RDFs in the coexistence region indicates 
that the typical values of the nearest-neighbour distance in all phases fall well within the 
repulsive shoulder of the pair potential, as expected for any highly-compressed system. 

An inspection of Table 2 shows that, as temperature grows, the difference P FS — P S s 
gradually reduces until it would cross zero at T* xs 19, as being suggested by a power-law 
extrapolation beyond T* = 16 of Pss(T) and Pps(T). Originally, it was just this evidence 
that led us to include in our calculations two other isotherms at T* = 20 and 25, in order to 
see whether a BCC-FCC phase transition eventually shows up at very high temperatures. 
Indeed, when T* is as high as 20, we find that, on increasing pressure, the sequence of 
stable phases of the exp-6 model is fluid-BCC-FCC, with transitions at (3Pfs = 83.517 r~ 3 
(fluid-to-BCC transition) and at (3P S s = 83.575 r m 3 (BCC-to-FCC transition). Actually, at 
this temperature the BCC window is so narrow that the fluid-BCC-FCC triple temperature 
should be very close to T* = 20. At T* = 25, the interval of stability of the BCC phase is 
much wider (f» 6 GPa). Table 3 reports thermodynamic data for T* = 20 and 25, whereas 
the difference in chemical potential between the various phases is shown in Fig. 4 for T* = 25 
as a function of pressure. Hence, the BCC phase becomes eventually stable, but at much 

14 



higher temperatures than estimated in Ref. [5]. 

In the end, the P-T phase diagram of the exp-6 model would appear as in Fig. 5, where 
temperature and pressure are reported in Xe units. In this picture, our fluid-solid coexistence 
data are contrasted with the coexistence loci of Belonoshko et al. and with the Simon-Glatzel 
equation of state for Xe, as deduced from that of Ne by the simple rescaling of pressure and 
temperature values proposed by Vos et al. [29]. Compared to the findings of Ref. [5], freezing 
of the exp-6 system occurs for slightly lower pressures in our simulation; more important, 
our fluid-FCC-BCC triple point lies very far from the location indicated by Belonoshko et 
al. 

To better understand the nature of the differences between our results and those of 
Belonoshko et al, we have plotted in Fig. 6 the exp-6 phase diagram on the T-V plane. It 
can be appreciated from this picture that small differences in the volumes of the coexisting 
fluid and FCC solid are found at all temperatures. On the contrary, isobaric paths are pretty 
much the same (see Fig. 11 of Ref. [5] for a comparison), indicating that where we deviate 
from Belonoshko et al. is essentially in the definition of a phase-coexistence condition, not 
in the calculation of basic statistical averages. 

Where we radically contrast with Belonoshko et al. is really in the position of the BCC- 
FCC transition line, with ours running much higher in temperature and pressure. To assess 
how much distant is our prediction from these authors, we note that, for e.g. T* = 16, 
they report a reduced pressure of about 70 GPa (P*/T* ~ 120) at the BCC-FCC transition 
while, at the same pressure, our PA/j, is about —0.040 (see Fig. 2), i.e., appreciably larger, 
in absolute terms, than the size of statistical errors. 

IV. DISCUSSION 

Having devoted the last Section to a plain presentation of our results, we now reconsider 
more critically their implications, especially in relation to the numerical precision of the 
simulation data. 
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The crucial quantity to look at is, obviously, the chemical potential. One source of error 
in its calculation follows from the finite size of the simulated system, a problem further 
complicated by the impossibility to compare FCC and BCC solids with equal numbers 
of particles. In point of principle, this would force us to some sort of extrapolation to 
N = oo which, however, is computationally demanding. Instead, we have decided to simulate 
fluid, FCC, BCC, and HCP samples of similar size, which are also sufficiently large that no 
significant finite-size error would be made in tracing the coexistence lines. Indeed, we have 
already demonstrated in the previous Section - see the comment to Fig. 2 - that this kind 
of error is not very important. 

Further errors in estimating the relative stability of two distinct solid phases are im- 
putable to the limited precision with which we compute the excess free energy of the refer- 
ence solid states in Table 1 as well as the pressure values along an isothermal path. Looking 
at Table 1, the overall error on the free-energy difference between BCC and FCC in units 
of /cgT is of about three units on the third decimal place, and comparable is the statistical 
error accompanying the values of f3P/p for each phase. The error associated with the Pfi 
gap between the fluid and a solid phase is again of 4 4- 5 x 1CT 3 , but the rate at which this 
quantity varies as a function of pressure is much larger (this implies that the location of 
the fluid-solid transition is much better defined numerically than is the solid-solid transition 
point). Summing up, we expect that the typical statistical error affecting the BCC-FCC 
PA.fi is of about 1CT 2 . If this is true, it is evident from the tabulated values of PAfi that we 
cannot definitely rule out the possibility that the BCC solid becomes stable at temperatures 
lower than T* = 20, and even as small as 11. Only if we turn to much longer Monte Carlo 
runs than hereby considered, we can hope to reduce drastically the width of the error bars. 
Anyway, even in the worst case, we can safely infer from our data that the BCC phase may 
only be stable in a narrow slice of few GPa adjacent to the fluid-BCC coexistence locus, a 
pressure interval much narrower than predicted by Belonoshko et al. 

The caution expressed above is probably overstated: it is a fact that all the curves 
plotted in Figs. 2 and 4 are very smooth, which suggests that the statistical noise underlying 
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the profile of /3A/J, is quite smaller than the maximum estimated above. In this case, the 
conclusions drawn in the previous Section just on the basis of the average behaviour of A/x 
are substantially correct, and the BCC phase would really be metastable below T* 20. 
If we believe this, a stable BCC phase would first appear in Xe at so high a temperature 
4500 K) that one can even wonder whether the BCC-FCC transition in Xe is preempted 
by quantum effects. 

As a matter of fact, our simulation results do not match the experimental data of Boehler 
et al. at the highest temperatures. Even assuming that the exp-6 system is a very good 
representation of Xe in the HT/HP regime, the problem could be in fact with the experiment, 
which might be far from realizing hydrostatic conditions. In fact, it is generally believed 
that, in a laser-heated DAC experiment, the stress state within the diamond cell may not 
be hydrostatic; moreover, the sample can even exhibit considerable shear stresses [30]. In 
consideration of the small difference in chemical potential between BCC, FCC, and HCP, it 
could then be that what is experimentally recognized as solid is in fact a mixture of BCC 
and FCC/HCP crystallites. 

V. CONCLUSIONS 

In this paper, we have reported on the results of an extensive NVT Monte Carlo simula- 
tion of the exp-6 model, an effective-interaction model that is thought to provide a realistic 
description of the thermal properties of rare gases under extreme, high-temperature/high- 
pressure conditions. We have plotted the "exact" phase diagram of this system upon com- 
bining the method of thermodynamic integration with fully-fledged free-energy calculations 
both in the fluid (by the Widom method) and in the solid phase (by the Frenkel-Ladd 
method). 

The aim of this effort was to point out some aspects concerning the uncertain status 
of the solid phase of Xenon at high densities. In a previous molecular-dynamics study of 
the exp-6 system, with parameters being appropriate to Xe, the two-phase method was 
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employed for simulating the thermodynamic coexistence between fluid and solid [5,10]. This 
study gave evidence of a stable BCC phase in a narrow range of pressures above 25 GPa, 
when temperature exceeds 2700 K, i.e., near to where a laser-heated, diamond- anvil-cell 
experiment finds a cusp on the freezing line of Xe. In fact, our free-energy calculations 
partially contradict such findings, showing that Xe freezes directly into a FCC solid, the 
BCC phase becoming stable only at much higher temperatures (above 4500 K), provided 
that the exp-6 modelization is still valid for Xe in these extreme conditions; moreover, the 
HCP solid is practically as stable as the FCC one. 

A clue to understand this disagreement is the very small difference in Gibbs free energy, 
for T > 2500 K, between the BCC and FCC solids of equal pressure and similar size. The 
little advantage of FCC over BCC at freezing could be the reason for the apparent stability 
of the BCC phase in the former numerical study of the exp-6 model. In our simulation study, 
we find that, at sufficiently high temperature, there is a narrow range of pressures where 
the BCC solid is more stable than the FCC, but less stable than the fluid. Only beginning 
from T* = 20, the BCC phase gains true thermodynamic stability in an interval of pressures. 
Obviously, this careful monitoring of the relative stability of the various phases as a function 
of both temperature and pressure would simply be impossible without the knowledge of their 
respective chemical potentials, and this is the ultimate reason for preferring exact free-energy 
calculations to other methods. 

As far as the experiment is concerned, it is likely that the cusp-like feature on the Xe 
freezing line is just an experimental artefact, which could be due to the increasing difficulty, 
as the freezing density progressively grows, in establishing hydrostatic conditions, a problem 
also worsened by the existence of two different solid phases that so closely compete with 
each other for stability. 
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TABLE CAPTIONS 



Table 1 : Excess free energy per particle in units of k B T for a number of exp-6 solid 
states. For each state and solid type, also shown within square brackets is the value of 
c* = cr 2 m je that intervenes the Frenkel-Ladd calculation: for the chosen c*, 5i?| in ap- 
proximately matches the average square deviation of an exp-6 particle from its position 
in the perfect crystal. 

Table 2 : Phase-transition data for T* = 4.25,8.15, 12.77, and 16. With the exception of 
T* = 4.25, simulation data refer to N = 1372 (fluid and FCC) and N = 1458 (BCC). 
For T* = 4.25, N = 864 for fluid and FCC, no BCC phase was considered. From 
left to right, values of (3Pss (BCC-to-FCC "virtual" transition, falling inside the fluid 
region of the phase diagram), /3Pps (fluid-to-FCC transition), (3fJ,ps (common value of 
(3 fx for the coexisting fluid and FCC solid), (3Alm fs (reduced chemical potential of FCC 
relative to BCC at P = Pfs), Pi (freezing density), and p m (FCC melting density) 
are shown. These quantities were arbitrarily round-off at the third decimal digit (the 
fourth digit only for /5A/ips)- However, the error accompanying them would usually 
be larger, originating from the limited precision of both the Monte Carlo data and the 
free-energy values in Table 1. 

Table 3 : Phase-transition data for T* = 20 and 25. From left to right, values of /3Pfs 
(fluid-to-BCC transition), PPss (BCC-to-FCC transition), pf (freezing density), p m 
(BCC melting density), Pbcc (BCC density at the BCC-FCC transition), and p F cc 
(FCC density at the BCC-FCC transition) are shown. 
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FIGURE CAPTIONS 



Fig. 1 : Frenkel-Ladd calculation of solid free energies: FCC solid at p* = 3.5 and T* = 8.15. 
Top: the integrand of Eq. (2.10) is plotted in the panel above for iV = 1372 (the 
continuous line is a spline interpolant of the data points and the error bars, which are 
also shown, are much smaller than the symbols size). In the panel below, the finite-size 
effect is demonstrated for (3 (AV) A /N through the difference between N = 500 and 
N = 1372 (A) and between N = 864 and N = 1372 (O)- Bottom: the values of 
Pfe*(N) + In N/N (for N = 500,864, and 1372) scale linearly with N" 1 for large N, 
as been conjectured in Ref. [26]. We have verified that the same type of scaling holds 
for all of the solid-state points in Table 1. 

Fig. 2 : Reduced chemical potential of the FCC solid relative to BCC and to HCP for 
T* = 16. Upon plotting the difference in (3fi between FCC and BCC (N = 1372 vs. 
N = 1458, continuous line; N = 864 vs. N = 1024, dashed line), we find that, as 
pressure grows, the BCC phase loses stability to the advantage of FCC. All the plotted 
lines are linear-spline interpolants of the data points. The dotted line is the quantity 
P(Hfcc ~ /Uhcp) (N = 1372 for the FCC solid; N = 1440 for the HCP solid). The two 
couples of vertical lines mark, from left to right, the position of the "virtual" transition 
between BCC and FCC and the fluid-FCC phase transition. 

Fig. 3 : Mechanical equation of state for T* = 16. The fluid branch (O and continuous 
line) and the FCC branch (□ and dotted line) are plotted for N = 1372 particles, 
while the BCC branch (A and dashed line) is for a system of 1458 particles. In the 
inset, we zoom into the transition region: the freezing and melting densities (signalled 
by vertical lines) are located where the horizontal line at (3P FS crosses the fluid and 
FCC equations of state. 

Fig. 4 : Reduced chemical potential of the FCC (N = 1372) phase relative to BCC (N = 
1458, continuous line) and HCP (N = 1440, dotted line), and of the fluid phase 
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(N = 1372) relative to BCC (N = 1458, dashed line), for T* = 25. Upon increasing 
pressure, we observe a fluid-BCC-FCC sequence of phases. 

Fig. 5 : HT/HP phase diagram of the exp-6 model (temperature and pressure are rescaled 
using Ross parameters for Xe). Besides our data (full and open dots), we plot the 
coexistence loci for the exp-6 model as been calculated by Belonoshko et al. (dotted 
lines), and the Xe freezing line as been derived - after suitable rescaling of temperature 
and pressure values - from the Simon-Glatzel fit of Ne data by Vos et al. (dashed line). 

Fig. 6 : Volume-temperature phase diagram of Xe as drawn from our numerical simulation 
of the exp-6 model. Present data for the melting and the freezing volumes (□ and 
continuous lines) and for the BCC-FCC coexistence volumes (A and straight lines) 
are compared with the estimates by Belonoshko et al. ( x and dotted lines, as got from 
Fig. 11 of Ref. [5]). Following Fig. 11 of the cited reference, we have also plotted as 
circles some isobars (open circles, fluid; full dots, FCC solid). 
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TABLE 1 



PJex 


r UL (iv = lot Z) 


U(~^T> (AT 1 A Al~\\ 

rivjr [1\ = 144U ) 


TiCC (AT 1 A KQ\ 


p* = 4,T* = 4.25 


38.879(1) [8500] 




39.446(2) [5100] 


p* = 3.5 , T* = 8.15 


16.258(1) [5400] 




16.408(2) [3000] 


p* = 4,T* = 12.77 


17.216(1) [7700] 




17.316(2) [5300] 


p* = 5 , T* = 16 


25.654(1) [15000] 


25.659(1) [14500] 


25.784(2) [9100] 


p* = 5.5, T* = 20 


27.289(1) [19600] 


27.295(1) [19400] 


27.397(2) [13600] 


p* = 5 , T* = 25 


18.489(1) [13800] 


18.494(1) [14300] 


18.532(2) [10100] 
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TABLE 2 
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P-TFS [T m ) 


PUfs 




Pf l r m J 


Pm m ) 


4.25 




28.897 


18.774 




1.980 


2.055 


8.15 


40.979 


46.883 


25.565 


-0.0124 


2.496 


2.571 


12.77 


57.445 


62.238 


29.519 


-0.0060 


2.946 


3.022 


16 


69.130 


71.877 


31.552 


-0.0029 


3.219 


3.295 
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TABLE 3 









Pf (r~ 3 ) 


Pm (r~ 3 ) 


Pbcc (r~ 3 ) 


Pfcc (r~ 3 ) 


20 


83.517 


83.575 
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3.611 
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96.901 
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4.048 
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